require(stats)
require(graphics)
#library("ggplot2")
require(ggplot2)
#library("grDevices")
require(grDevices)
#library("psych")
require(psych)
library(gridExtra)
library(grid)

# Please make sure to change the path to where the appropriate dataset is located on your computer

sc.val <-  read.csv("PNAS_VAL.csv") 


# Columns

sponsor.require <- 1
p.tourist.visa <- 2
p.fakedocs <-  3
k <- 4
my.seed <- 5

legaltourist.attempts <- c(6:25)
ability.legal.cols <- c(26:86)
ability.fakedocs.cols <- c(87:147)
ability.work.cols <- c(148:208)

tot.mig <- c(209:269)
legal.mig <- c(270:330)
fd.mig <- c(331:391)
ut.mig <- c(392:452)
stud.mig <- c(453:513)
hs.mig <- c(514:574)
ls.mig <- c(575:635)
fam.mig <- c(636:696)

legal.start <- 232
ut.start <- 9
fd.start <-  9
tot.start <-  250

####

table(sc.val[,k])

sc.val$tot.legal <- as.vector(sc.val[,legal.mig[1:12]] + legal.start)
sc.val$tot.ut <- as.vector(sc.val[,ut.mig[1:12]] + ut.start)
sc.val$tot.fd <- as.vector(sc.val[,fd.mig[1:12]] + fd.start)

sc.val$tot.illegal <- sc.val$tot.ut + sc.val$tot.fd

sc.val$all.mig <- sc.val$tot.legal + sc.val$tot.ut + sc.val$tot.fd

sc.val$prop.leg <- sc.val$tot.legal / sc.val$all.mig

sc.val$prop.ut <- sc.val$tot.ut / sc.val$all.mig
sc.val$prop.fd <- sc.val$tot.fd / sc.val$all.mig
sc.val$prop.ill <- sc.val$prop.ut + sc.val$prop.fd

percentile.ill.k1 = quantile(sc.val$prop.ill[which(sc.val[,k] == 1),12], c(0.0275, 0.975))
percentile.leg.k1 = quantile(sc.val$prop.leg[which(sc.val[,k] == 1),12], c(0.0275, 0.975))
percentile.ill.k2 = quantile(sc.val$prop.ill[which(sc.val[,k] == 2),12], c(0.0275, 0.975))
percentile.leg.k2 = quantile(sc.val$prop.leg[which(sc.val[,k] == 2),12], c(0.0275, 0.975))
percentile.ill.k3 = quantile(sc.val$prop.ill[which(sc.val[,k] == 3),12], c(0.0275, 0.975))
percentile.leg.k3 = quantile(sc.val$prop.leg[which(sc.val[,k] == 3),12], c(0.0275, 0.975))


tab.leg.illeg.k1 <- cbind(sc.val$prop.leg[which(sc.val[,k] == 1),12], sc.val$prop.ill[which(sc.val[,k] == 1),12])
tab.leg.illeg.k2 <- cbind(sc.val$prop.leg[which(sc.val[,k] == 2),12], sc.val$prop.ill[which(sc.val[,k] == 2),12])
tab.leg.illeg.k3 <- cbind(sc.val$prop.leg[which(sc.val[,k] == 3),12], sc.val$prop.ill[which(sc.val[,k] == 3),12])

tab.census <- c(0.85, 0.15)

sim.leg.illeg.k1 <- colMeans(tab.leg.illeg.k1)
sim.leg.illeg.k2 <- colMeans(tab.leg.illeg.k2)
sim.leg.illeg.k3 <- colMeans(tab.leg.illeg.k3)

legileg.k1 <- rbind(tab.census,sim.leg.illeg.k1)
legileg.k2 <- rbind(tab.census,sim.leg.illeg.k2)
legileg.k3 <- rbind(tab.census,sim.leg.illeg.k3)

legileg.k <- rbind(tab.census,sim.leg.illeg.k1,sim.leg.illeg.k2,sim.leg.illeg.k3)


sc.val$prop.mig.k1 <-  sc.val$all.mig[which(sc.val[,k] == 1),12]/1166
sc.val$prop.mig.k2 <-  sc.val$all.mig[which(sc.val[,k] == 2),12]/1166
sc.val$prop.mig.k3 <-  sc.val$all.mig[which(sc.val[,k] == 3),12]/1166

percentile.mig.k1 = quantile(sc.val$prop.mig.k1, c(0.0275, 0.975))
percentile.mig.k2 = quantile(sc.val$prop.mig.k2, c(0.0275, 0.975))
percentile.mig.k3 = quantile(sc.val$prop.mig.k3, c(0.0275, 0.975))

tab.prop.mig.k <- c(0.34, mean(sc.val$prop.mig.k1), mean(sc.val$prop.mig.k2), mean(sc.val$prop.mig.k3))

k.names.2 <- c("Census", "ABM, k=1","ABM, k=2", "ABM, k=3","Census", "ABM, k=1","ABM, k=2", "ABM, k=3")
k.names <- c("Census", "ABM, k=1","ABM, k=2", "ABM, k=3")
k.colors <- c("gray42",rgb(0,1,0.5, alpha = 0.9),rgb(0,1,0.5, alpha = 0.6), rgb(0,1,0.5, alpha = 0.4))


pdf("FigC4", height=5, width=5)
nf <- layout(matrix(c(1,2), 1, 2, byrow = TRUE), widths=c(1,1))

# bar graph
par(pty="s", mar=c(6,4,3,2), xpd=TRUE)
barplot(legileg.k, col= k.colors, ylim = c(0,1), horiz = F, beside = T, ylab = "", 
        xlab = "", axes = T, border = F, names.arg = k.names.2, cex.names = 0.8, las = 3,cex.axis = 1)
segments(2.5,percentile.leg.k1[1],2.5,percentile.leg.k1[2], lty=1)
segments(7.5,percentile.ill.k1[1],7.5,percentile.ill.k1[2], lty=1)
segments(3.5,percentile.leg.k2[1],3.5,percentile.leg.k2[2], lty=1)
segments(8.5,percentile.ill.k2[1],8.5,percentile.ill.k2[2], lty=1)
segments(4.5,percentile.leg.k3[1],4.5,percentile.leg.k3[2], lty=1)
segments(9.5,percentile.ill.k3[1],9.5,percentile.ill.k3[2], lty=1)

mtext("(a) Legal       Illegal", side=3, line=1, cex=1, font = 4)
mtext("Proportion of Migrants", side=2, line=2, cex=1)

### Volume / Population of Jamaica
par(pty="s", mar=c(6,4,3,2), xpd=TRUE)
barplot(tab.prop.mig.k, col=k.colors, beside = TRUE, space = 0, border = F, ylim=c(0,1),
        cex.axis = 1, names.arg = k.names, cex.names = 0.8, las = 3,
        ylab = "", 
        xlab = "" )
segments(1.5,percentile.mig.k1[1],1.5,percentile.mig.k1[2], lty=1)
segments(2.5,percentile.mig.k2[1],2.5,percentile.mig.k2[2], lty=1)
segments(3.5,percentile.mig.k3[1],3.5,percentile.mig.k3[2], lty=1)

mtext("Migrants : Pop. Jamaica", side=2, line=2.5,cex=1)
mtext("(b) Volume", side=3, line=1, cex=1, font = 4)
dev.off()
